/*
 *  @file    graph.h
 *  @author  Yi He
 *  @brief   template for graph and graph-related algorithms.
 *  @version $Id: LICENSE.h 1615 2010-10-20 19:16:35Z sobolevnrm $
 *  @attention
 *  @verbatim
 * This file is part of Jamming
 *
 * Copyright (C) 2012 - Yi He
 *
 * Jamming is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Jamming is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Jamming. If not, see <http://www.gnu.org/licenses/>.
 * @endverbatim
 */


#ifndef __GRAPH_H__
#define __GRAPH_H__
#include<map>
#include<vector>
#include<utility>
#include<iostream>
#include<queue>
#include<exception>
#include<limits>
#include<stack>
#include "debug.h"
#include "tree/disjoint_sets.h"

namespace graph
{
/***
*  Graph data structure using adjacent list.
*/
template<class T>
class adjacent_graph  {
        std::vector<std::map<int,T> > g;
        public:
        void addedge(int i,int j,T v) {
            if(i>=g.size()) g.resize(i+1,std::map<int,T>());
            //TODO allow i==j? 
            if(v!=0) g[i][j]=v;
        }
        inline int numvertex() {
            return g.size();
        }
        inline T& weight(int i,int j) {
            //TODO currently, visit noexist (i,j) will generate a element with value 0(T()). It is better to create a weight function that generate no empty space in graph.
              return g[i][j];
        }
        inline T& operator() (int i,int j) {
            if(i < g.size()) return g[i][j];
            else {
                g.resize(i+1);
                return g[i][j];
            }
        }
        //input&output
        //import from a array/matrix with N lines and M columns.
        void import(T* arr,int N,int M) {
            for(int i=0;i<N;i++)
                for(int j=0;j<M;j++) {
                    addedge(i,j,*arr);
                    arr++;
                }
        }
        //read from a stream
        template <class inputstream>
        void read(inputstream& in,int N,int M) {
            T v;
            for(int i=0;i<N;i++){
               std::string l;
              for(int j=0;j<M;j++) {
                in>>v;
                addedge(i,j,v);
              }
              getline(in,l);
            }
        }
        template <class outputstream>
        void print_edges(outputstream &out) {
            int w=g.size();
            for(int i=0;i<w;i++) {
                out<<i<<" : ";
                for(typename std::map<int,T>::iterator it=g[i].begin();it!=g[i].end();it++) {
                    if(it->second != T()) out<<"("<<it->first<<","<<it->second<<")";
                    else g[i].erase(it);
                }
                out<<std::endl;
            }
        }
        template <class outputstream>
        void print_map(outputstream &out) {
            int w=g.size();
            for(int i=0;i<w;i++) {
              for(int j=0;j<w;j++)
                out<<g[i][j]<<" ";
               out<<std::endl;
            }
        }
        //iterator 
        class edge_iterator {
            adjacent_graph<T>* pG; // graph pointed to.
            typename std::map<int,T>::iterator pC; //current position
            int pV; //source vertex pointed to. 
            //condense zero edges
            inline void condense() {
                    //Remove edge that equals T() (-> 0 -> no connection).
                    while(pC->second==T()) {
                        if(pC==pG->g[pV].end()) break;
                        pG->g[pV].erase(pC++);
                    }
            }
            public:
                edge_iterator();
                edge_iterator(adjacent_graph<T>* p,int v,typename std::map<int,T>::iterator c) : pG(p),pV(v),pC(c) {}
                inline int vertex() {
                    condense();
                    return pC->first;
                }
                inline T length() {
                    condense();
                    return pC->second;
                }
               edge_iterator& operator++ () {
                    if(pC->second!=T()) pC++;
                    else do {
                    //condense zero edges
                        if(pC==pG->g[pV].end()) return *this;
                        pG->g[pV].erase(pC++);
                    } while(pC->second==T());
                    return *this;
                }
                //TODO implement it++? (force to use ++it so that there is no additional performance cost?)
                bool operator!= (const edge_iterator& ref) {
                    condense();
                    if(pG==ref.pG && pV==ref.pV && pC==ref.pC) return false;
                    else return true;
                }
                bool operator==(const edge_iterator& ref) {
                    condense();
                    if(pG==ref.pG && pV==ref.pV && pC==ref.pC) return true;
                    else return false;
                }
                //TODO implement operator* & operator->
                std::pair<const int,T> & operator* () {
                    condense();
                    return *pC;
                }
                typename std::map<int,T>::iterator operator->() {
                    condense();
                    return pC;
                }
        };
        edge_iterator first_edge(int v) {
            return edge_iterator(this,v,g[v].begin());
        }
        edge_iterator end_edge(int v) {
            return edge_iterator(this,v,g[v].end());
        }
}; 


/*!
* Graph explorer
*/
// exceptions for explorer
class GraphHasLoop : public std::exception {};
// Helper classes for explorer
//Conditions to end dfs/bfs search 
template <class T>
class reached_at {
        int t;
        public:
            reached_at(int ter) : t(ter) {};
        bool operator() (int v,int prev, bool hasOutEdge) {
                return (v==t);
            }
        bool operator() (std::vector<T>& dist) {
            if(dist[t]<std::numeric_limits<T>::max()) return true;
            return false;
        }
};
bool visited_all(int v,int prev, bool hasOutEdge) { return false;}
class terminate_at{
        int& t;
        public:
            terminate_at(int& ter) : t(ter) {};
        bool operator() (int v,int prev, bool hasOutEdge) {
                t=v;
                return !hasOutEdge;
            }
};

//conditions to visit an edge in dfs/bfs search 
template <class T>
class hasedge {
    public:
    bool operator()(int v,int prev,T& edge) { return (edge!=T());}
};


template <class T,class container=class adjacent_graph<T> >
class explorer {
    T Inf;
    std::vector<std::vector<T> > dist;
    std::vector<std::vector<int> > prev;
    container& g; 
    public:
    explorer() {
        Inf=std::numeric_limits<T>::max();
    }
    explorer(container &c) : g(c) {
        Inf=std::numeric_limits<T>::max();
    }
    /*! 
     * bfs walker 
    */
    template <class _Vertex_found_func ,class _Edge_check_func >
    bool bfs_until(int s,_Vertex_found_func found,_Edge_check_func edge_check) {
        int u,v;
        std::queue<int > q; 
        if(s>=dist.size()) dist.resize(s+1);
        if(s>=prev.size()) prev.resize(s+1);
        prev[s]=typename std::vector<int> (g.numvertex(),-1);
        typename std::vector<bool> visited(g.numvertex(),0);
        visited[s]=true;
        q.push(s);
        while(!q.empty()) {
            u=q.front();
            q.pop();
            bool hasOutEdge=false;
            for(typename container::edge_iterator it=g.first_edge(u);it!=g.end_edge(u);++it) {
                v=it.vertex();
                T w=it.length(); //TODO edge requirment-> comp() function.
                if(!visited[v]&&edge_check(v,u,w)) {
                    prev[s][v]=u;visited[v]=true;
                    hasOutEdge=true;
                    q.push(v);
                }
            }
            if(found(u,prev[s][u],hasOutEdge)&&u!=s) return true;
        }
        return false;
    }
    template <class _Vertex_found_func >
    bool bfs_until(int s,const _Vertex_found_func& found) {
        return bfs_until(s,found,hasedge<T>());
    }
    bool bfs(int s,int t=-1) {
        return bfs_until(s,reached_at<T>(t),hasedge<T>());
    }
    // dfs walker 
    template <class _Vertex_found_func ,class _Edge_check_func >
    bool dfs_until(int s,_Vertex_found_func found,_Edge_check_func edge_check) {
        if(s>=dist.size()) dist.resize(s+1);
        if(s>=prev.size()) prev.resize(s+1);
        std::stack<int> q;
        prev[s]=std::vector<int> (g.numvertex(),-1);
        std::vector<bool> visited(g.numvertex(),0);
        q.push(s);
        visited[s]=true;
        while(!q.empty()) {
            int u=q.top(); 
            bool hasOutEdge=false;
            for(typename container::edge_iterator it=g.first_edge(u);it!=g.end_edge(u);++it) {
                int v=it.vertex();
                T w=it.length(); //TODO edge requirment-> comp() function.
                if(!visited[v]&&edge_check(v,u,w)) {
                    prev[s][v]=u;visited[v]=true;
                    q.push(v);
                    hasOutEdge=true;
                }
            }
            if(found(u,prev[s][u],hasOutEdge)&&u!=s) return true;
            if(!hasOutEdge) q.pop();
        }
        return false;
    }
    template <class _Vertex_found_func >
    bool dfs_until(int s,const _Vertex_found_func& found) {
        return dfs_until(s,found,hasedge<T>());
    }
    bool dfs(int s,int t=-1) {
        return dfs_until(s,reached_at<T>(t),hasedge<T>());
    }
    //shortest distance walkers
    /**
    * dijkstra algorithm, find the shortest distance between source and destiny.
    * 
    */
    template <class _Vertex_found_func >
    bool dijkstra_until(int s,_Vertex_found_func found) {
        T w;
        if(s>=dist.size()) dist.resize(s+1);
        if(s>=prev.size()) prev.resize(s+1);
        dist[s]=std::vector<T>(g.numvertex(),Inf);
        prev[s]=std::vector<int>(g.numvertex(),-1);
        //if(s>=g.numvertex()) throw GraphHasLoop();
        dist[s][s]=T();
        typename std::priority_queue<std::pair<T,int> ,std::vector<std::pair<T,int> >, std::greater<std::pair<T,int> > > q;
        q.push(std::pair<T,int>(T(),s));
        std::vector<bool> visited(g.numvertex(),false);
        while(!q.empty()) {
            std::pair<T,int> tV=q.top(); q.pop(); //top vertex.
            int u=tV.second;
            visited[u]=1;
            for(typename container::edge_iterator it=g.first_edge(u);it!=g.end_edge(u);++it) {
                int v=it.vertex(); w=it.length(); 
                //TODO ? what about formula description? redefine its length?
                if(visited[v]) continue;
                //relax
                if( dist[s][v] > dist[s][u]+w) {
                     dist[s][v]= dist[s][u]+w; prev[s][v]=u;
                }
               q.push(std::pair<T,int>(dist[s][v],v));
            }
        }
        return false;
    }
    bool dijkstra(int s,int t=-1) {
        return dijkstra_until(s,reached_at<T>(t));
    }
    /**
    * SPFA algorithm.  A modification of bellman algorithm. Allows negative edge.
    *
    */
    template <class _Vertex_found_func >
    bool spfa_until(int s,_Vertex_found_func found) {
        if(s>=dist.size()) dist.resize(s+1);
        if(s>=prev.size()) prev.resize(s+1);
        dist[s]=std::vector<T>(g.numvertex(),Inf);
        prev[s]=std::vector<int>(g.numvertex(),-1);
        dist[s][s]=0;
        std::deque<int> q;
        q.push_back(s);
        std::vector<bool> inQ(g.numvertex(),0);
        while( q.size()!=0) {
            int u=q.front(); q.pop_front(); inQ[u]=false;
            //if(u==t) return true;
            //TODO negative loop check?
            for(typename container::edge_iterator it=g.first_edge(u);it!=g.end_edge(u);++it) {
                int v=it.vertex(); T w=it.length(); //TODO what about formula description?
                //relax
                if( dist[s][v] > dist[s][u]+w) {
                    dist[s][v]=dist[s][u]+w; prev[s][v]=u;
                     if(!inQ[v]) {
                        inQ[v]=true;
                        if(q.size()!=0) {
                            if(dist[s][v]>dist[s][q.front()]) q.push_back(v);
                            else q.push_front(v);
                        } else {
                            q.push_back(v);
                        }
                     }
                }
            }
        }
        if(found(dist[s])) return true;
        return false;
    }
    bool spfa(int s,int t=-1) {
        return spfa_until(s,reached_at<T>(t));
    }
    /**
    * flyod algorithm, find the shortest distance between every vertex pairs.
    *
    */
    void flyod() {
        dist=std::vector<std::vector<T> > (g.numvertex(), std::vector<T>(g.numvertex(),Inf) );
        //prev=std::vector<std::vector<int> > (g.numvertex(), std::vector<int>(g.numvertex(),-1) );
        for(int i=0;i<g.numvertex();i++) 
          for(int j=0;j<g.numvertex();j++) {
            if(i==j) dist[i][j]=0;
            else if(g.weight(i,j)!=0) dist[i][j]=g.weight(i,j);
            else dist[i][j]=Inf;
        }
        for(int k=0;k<g.numvertex();k++) 
           for(int i=0;i<g.numvertex();i++) 
            if(i!=k) 
              for(int j=0;j<g.numvertex();j++) {
                if(i!=j&&j!=k && dist[i][k]< Inf && dist[k][j]< Inf \
                    && dist[i][k]+dist[k][j]<dist[i][j]) {
                    dist[i][j]=dist[i][k]+dist[k][j];
                    //path[i][j]=k;
                }
              }
    }
    //output distance
    inline T distance(int s,int t) {
        return dist[s][t];
    }
    // empty the space occupied by explorer
    void clear() {
        dist.resize(0);
        prev.resize(0);
    }
    //path iterator
    class path_iterator {
        container& g;
        std::vector<int>& prev;
        int cur_pos; //point to current position;
        public:
            path_iterator(container& g,std::vector<int>& chain,int pos) :g(g),prev(chain),cur_pos(pos) {};
            bool operator!=(const path_iterator& ref) {
                if(cur_pos==ref.cur_pos) return false;
                else return true;
            }
            path_iterator& operator --() {
                cur_pos=prev[cur_pos];
                return *this;
            }
            int vertex() {
                return cur_pos;
            }
            int operator*() {
                return cur_pos;
            }
            T& operator-() {
                return g.weight(prev[cur_pos],cur_pos);
            }
            T& operator+() {
                return g.weight(cur_pos,prev[cur_pos]);
            }
    };
    path_iterator path_begin(int s,int t=-1) {
        return path_iterator(g,prev[s],s);
    };
    path_iterator path_end(int s,int t=-1) {
        if(t<0);
        return path_iterator(g,prev[s],t);
    };
};

template <class T,class container=class adjacent_graph<T> >
class connect {
    T Inf;
    container& g; //original graph.
    public:
    connect();
    connect(container& g) :g(g) {};
    //minimal span tree
    /**
    *  Kruskal algorithm. Used to find minimal span tree. Best case is a single tree, otherwise it will return a forest.
    */
    typedef std::pair<T,std::pair<int,int> > dist_vertex_pair;
    typedef std::multimap<T,std::vector<int> > span_tree;
     span_tree kruskal() {
        disjoint_sets<int> forest;
        typename std::priority_queue<dist_vertex_pair,std::vector<dist_vertex_pair>, std::greater<dist_vertex_pair> > q;
        //initialize 
        for(int i=0;i<g.numvertex();i++) {
            for(typename container::edge_iterator it=g.first_edge(i);it!=g.end_edge(i);++it) {
                q.push(dist_vertex_pair(it.length(),std::pair<int,int>(i,it.vertex())));
            }
            forest.addset(i);
        }
        int k=0,p1,p2;
        std::vector<T> sum(g.numvertex(),0);
        while(!q.empty()) {
            dist_vertex_pair t=q.top();q.pop();
            int u=t.second.first,v=t.second.second;
            if( (p1=forest.find(u))!=(p2=forest.find(v))) {
                sum[forest.merge(u,v)]=sum[p1]+sum[p2]+g.weight(u,v);
                k++;
            }
        }
        //organize for output 
        span_tree mst;
        disjoint_sets<int>::set_map r=forest.organize();
        for(disjoint_sets<int>::set_map::iterator it=r.begin();it!=r.end();it++) {
            if(it->second.size()>0) {
                mst.insert(typename std::pair<T,std::vector<int> > (sum[it->first],it->second));
            }
        }
        return span_tree(mst);
    }
    //Strong connected component
    /*
    *  tarjan algorithm  O(N+M)
    */
    class tarjan_visit {
        std::vector<int>& DFN;
        std::vector<int>& Low;
        std::stack<int> s;
        public:
        tarjan_visit(std::vector<int>& DFN,std::vector<int>& Low): DFN(DFN),Low(Low) {
//            DFN=std::vector<int> (g.numvertex(),std::numeric_limits<int>::max());
//            Low=std::vector<int> (g.numvertex(),std::numeric_limits<int>::max());
        };
        bool operator() (int v,int u, bool hasOutEdge) {
            std::cout<<"Visited "<<v<<std::endl;
            return false;
        }
    };
    void tarjan(int u) {
        //stack<int> s;
        //vector<bool> visited(g.numvertex(),0);
        std::vector<int> DFN(g.numvertex(),std::numeric_limits<int>::max());
        std::vector<int> Low(g.numvertex(),std::numeric_limits<int>::max());
	    DFN[u]=Low[u]=0;
	    explorer<T,container> walk(g);
        walk.dfs_until(u,tarjan_visit(DFN,Low));
    }
};

//maximum flow through the network
template <class T,class container=class adjacent_graph<T> >
class flow {
    T Inf;
    container& g; //original graph.
    public:
    flow() {
        Inf=std::numeric_limits<T>::max();
    }
    flow(container& c) : g(c) {
        Inf=std::numeric_limits<T>::max();
    }
    /**
    * EK algorithm for maximum flow. 
    *  O(V*E)
    */
    T EdmondsKarp(int s,int t) {
        T flow=0,d; //maximum flow & delta increment.
        int u;
        container* pRemain=new container;
        *pRemain=g;
        explorer<T,container> walk(*pRemain);
        while(walk.bfs(s,t)) {
           d=Inf;
           for(typename explorer<T,container>::path_iterator it=walk.path_end(s,t);it!=walk.path_begin(s,t);--it) d=std::min(d,-it);
           for(typename explorer<T,container>::path_iterator it=walk.path_end(s,t);it!=walk.path_begin(s,t);--it) {
                (-it)-=d;
                (+it)+=d;
            }
       flow+=d;
     }
      delete pRemain;
      return flow;
    }
    /**
    * dinic algorithm for maximum flow. 
    *  O(V^2 * E)
    */
    class dinic_build_level_til {
            std::vector<int>& level;
            int t;
        public:
            dinic_build_level_til(std::vector<int>& lvl,int end) :level(lvl),t(end) {}
            bool operator() (int v,int prev, bool hasOutEdge) {
                if(prev>=0) level[v]=level[prev]+1;
                else level[v]=0;
                if(v==t) return true;
                else return false;
            }
    };
    template <class _edge_type>
    class check_dinic_level {
            std::vector<int>& level;
    public:
        check_dinic_level(std::vector<int>& lvl) : level(lvl) {}
        bool operator()(int v,int prev,_edge_type& edge) { return (edge>T()&&level[v]==level[prev]+1);}
    };
    T dinic(int s,int t) {
        T flow=T(),d; //maximum flow & delta increment.
        int u;
        std::vector<int> level(g.numvertex(),0);
        container* pRemain=new container;
        *pRemain=g;
        explorer<T,container> walk(*pRemain);
        while(walk.bfs_until(s,dinic_build_level_til(level,t))) {
            d=Inf;
            while(walk.dfs_until(s,reached_at<int>(t),check_dinic_level<int>(level)) ) {
                    for(typename explorer<T,container>::path_iterator it=walk.path_end(s,t);it!=walk.path_begin(s,t);--it) d=std::min(d,-it);
                    for(typename explorer<T,container>::path_iterator it=walk.path_end(s,t);it!=walk.path_begin(s,t);--it) {
                        (-it)-=d;(+it)+=d;
                    }
                    flow+=d;
                }
        }
        delete pRemain;
        return flow;
    }
    /**
    * SAP(shortest-augment-path) algorithm for maximum flow
    */
    
    /**
    * Pre-flow algorithm for maximum flow
    */
    
};

template <class T,class container=class adjacent_graph<T> >
class bipartite {
    T Inf;
    container& g; //original graph.
    public:
    bipartite();
    bipartite(container& g) :g(g) {};
    /*
    * Hugarian Algorithm
    */
    class hg_matched {
        std::vector<int>& match;
        public:
            hg_matched(std::vector<int>& match) :  match(match) {}
            bool operator() (int v,int prev,bool hasOutEdge) {
                if(match[prev]<0) {
                    match[prev]=v;
                    return true;
                }
                return false;
            }
    };
    int hungarian() {
        int num_match=0;
        int u,t;
        std::vector<int> match(g.numvertex(),-1);
        explorer<T,container> walk(g);
        for(int s=0;s<g.numvertex();s++) 
            if (walk.bfs_until(s,hg_matched(match))) {
                num_match++;
            }
        return num_match;
    }
    /*
    *  Hopcroft-Karp Algorithm O(sqrt(n)*e)
    */
    int Hopcroft_Karp() {
    };
    /*
    *  Kuhn_Munkras Algorithm
    */
    int Kuhn_Munkras(){
    };
};

//cost Flow 
//template <class T,class container=class adjacent_graph<T> >
//class flow : {
//        //minCost Maxinum flow
//    //each edge has two value: cost & cap.
////    T EKminCostflow(int s,int t) {
////        T flow=0,d; //maximum flow & delta increment.
////        int u;
////        explorer<T,container> walk(*pG);
////        while(walk.spfa(s,t)) {
////           d=Inf;
////           for(u=t;u!=s&&u>=0;u=prev[u]) d=std::min(d,g.cap(prev[u],u));
////           for(u=t;u!=s&&u>=0;u=prev[u]) {
////                g.cap(prev[u],u)-=d;
////                g.cap(u,prev[u])+=d;
////                flow+=d*g.cost(u,prev[u]);
////            }
////     }
////      return flow;
////    }
//};
}; /* graph */


#endif /* __GRAPH_H__ */


